Exotic Ising dynamics in a Bose-Hubbard model 
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We explore the dynamical properties of a one-dimensional Bose-Hubbard model, where two different bosonic 
species interact via Feshbach resonance. We focus on the region in the phase diagram which is described by 
an effective, low-energy ferromagnetic Ising model in both transverse and longitudinal fields. In this regime, 
we numerically calculate the dynamical structure factor of the Bose-Hubbard model using the Time-Evolving 
Block Decimation method. In the ferromagnetic phase, we observe both the continuum of excitations and the 
bound states in the presence of a longitudinal field. Near the Ising critical point, we observe the celebrated E 8 
mass spectrum in the excited states. We also point out possible measurements which could be used to detect 
these excitations in an optical lattice experiment. 
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The search for emergent excitations arising from strong cor- 
relations has been very fruitful over the past few decades, 
from fractional quasiparticles in the fractional quantum Hall 
effect |H, to effective magnetic monopoles in spin ice Q. The 
one-dimensional (ID) transverse-field Ising model is another 
famous case showing collective excitations 0. It is a paradig- 
matic model for quantum phase transitions, hosting a critical 
point with central charge c=l/2 |4|. Zamolodichkov showed 
that, by perturbing such a critical theory with a Z 2 symmetry- 
breaking field, eight massive particles emerge in the excitation 
spectrum (5). These particles are the hallmark of an underly- 
ing E$ continuous symmetry, a very complex symmetry group 
discovered in mathematics, which attracts much interest in a 
wide community, see e.g. Ref. |6|. The ratio between the 
masses of the two lightest particles predicted by E% symme- 
try has been recently observed by Coldea et al. in neutron- 
scattering studies of the quantum magnet CoNb206 Q. 

The quest for complex many-body phenomena and uncon- 
ventional excitations has greatly benefited from advances in 
ultra-cold atomic gases 0[lO]|. By confining the atomic cloud 
to an optical lattice, it is now possible to explore strongly- 
interacting lattice models with an unprecedented degree of 
control and tunability (9). Important milestones include the 
observation of the superfluid to Mott insulator transition in 
bosonic systems ifTTIl and the antiferromagnetic ID Ising tran- 
sition on a tilted optical lattice \12\. Due to their long co- 
herence times and the possibility of tuning the parameters of 
the system, even as function of time, cold atoms in optical 
lattices also allow the study of non-equilibrium dynamics in 
these models [Q31 []31 - which is usually very difficult in a 
condensed-matter setting. A promising research area is the 
use of multi-component atomic mixtures as a route to effec- 
tive magnetic models, see e.g. Refs. fl"5lfT7ll . More specifi- 
cally, recent theoretical work has proposed mixtures of atoms 
and molecules near a Feshbach resonance [18 ] in a Mott state 
as a route to the ID Ising model fT9ti2T1l . Here, the Feshbach 
resonant coupling between different species plays the role of 
a tunable handle on quantum fluctuations. 

In this paper, we study the low-energy dynamical properties 
of a ID Bose-Hubbard model describing two different bosonic 



species coupled by Feshbach resonance in the Mott insulating 
regime. We find these to be captured by an appropriate dy- 
namical structure factor associated with the low-energy effec- 
tive Ising model. The defining characteristics of the broken- 
symmetry and disordered phases are clearly observed, both in 
the presence and absence of an effective longitudinal field. By 
tuning the parameters of the system close to the Ising quantum 
critical point, the excitation spectrum reveals the signatures of 
the emerging E$ symmetry. Furthermore, we point out pos- 
sible directions to detect these excitations in an optical lattice 
experiment. 

We consider the following pairing Bose-Hubbard Hamilto- 
nian, previously introduced in fT9U27lL 

U =J2 e * n ^ ~ X^^lcA+ia +H.c.) 

iot ia 

+ Yl ^Y~ n ^(n ia ' ~ <W) + gY,( b lm b ia b ia + H.C.), 

iota,' i 

(1) 

describing two species of bosons b ia on a ID lattice, where 
nia = b\ a b ia . Atoms are labeled by a = a while molecules 
are labeled by a = m. Here e a are on-site potentials, t a 
are hopping parameters between nearest-neighbour sites, and 
U aa t are on-site interactions. Two atoms form a molecule 
via 5-wave pairing, driven by Feshbach coupling g. The 
Feshbach interaction breaks the independent conservation of 
the number of atoms and molecules, but the total number 
N T = Xli( n m + 2n^ m ) is conserved. We work in the canon- 
ical ensemble, by keeping the total density pt=Nt/L fixed. 

The low-energy behaviour of the Hamiltonian Eq. ([T]) in the 
Mott regime with p T =2 is captured by an effective ID quan- 
tum Ising model fT9U2T1l . The effective spin degrees of free- 
dom are |ff) = |1; 0) and |J|) = |0; 2) in the basis of occupa- 
tion numbers \n a ; n m ), see Fig.[TJa). We truncate the Hilbert 
space to a maximum of three atoms and one molecule per site, 
which is already a good approximation to the canonical soft- 
core limit 1 21 1. This choice allows the hopping of atoms, even 
if a pair is already present on a site. The effective Ising model 
(up to an additive constant) is obtained via a strong-coupling 
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expansion around the small-hopping limit 

H ~ - J S i S Ui + hY,S!+Tj2sf + 0(f). (2) 



The spin operators have a direct physical interpreta- 
tion: Sf = (n irn — n ia /2)/2 = An^/2 measures the im- 
balance in the density of bosons at site i, while 
Sf = l/(2V2)[bl m b ia b ia + b im b\ a bl a ] accounts for fluctua- 
tions between the two species. At a qualitative level, the 
Ising exchange interaction J captures the motion of atoms and 
molecules, the longitudinal field h tunes an overall imbalance 
between the two species, and the transverse field T controls 
the fluctuations between the two species, cf. Fig.[TJa). 

We find the ground state of the full Hamiltonian Eq. ([I]) 
using a variant of the infinite density-matrix renormalization 
group (iDMRG) method (28| [29), which yields a matrix- 
product state (MPS) representation of the ground- state wave- 
function directly in the thermodynamic limit. We then use 
the infinite time-evolving block decimation (iTEBD) method 
l30l [311 to compute the transverse dynamical factor structure 
function S y (k,uj) of the effective spin representation at zero 
temperature. This is obtained by a Fourier transformation of 
the dynamical two-point correlation function 

cy(i,t) = (ip \si(t)s%(o)\<Po), O) 

found by evolving the infinite-size ground-state wave func- 
tion |^o) after applying the operator S y at a given site (32]- 
1351 . The sampled time is doubled by extrapolating C y {i,t) 
using linear prediction [36]. We stop the simulation once the 
"light-cone" of local correlations gets close to the boundary of 
a fixed observation window size (typically L w 200). Since 
the low-energy dynamics are set by the effective model, the 
"light-cone" and the entanglement entropy grow very slowly 
with time. Hence, we are able to reach extremely large times, 
up to tmax ~ 10 5 , while keeping the truncation error down 
to < 1 x 10 ~ 7 by only using an MPS for the ground state 
with bond dimensions xgs % 30. We have checked the con- 
vergence of our results by examining the dependence of the 
measured observables on \ and the Trotter time step At, set- 
tling on Xt = 60 and At = 0.1 — 0.2 for the simulations 
presented here. 

The ground-state properties of the pairing Bose-Hubbard 
model have been studied recently iflWTI l24l l25l While, 
previously, emphasis was placed on the phase diagram where 
the effective model is the antiferromagnetic Ising chain, here 
we focus on parameters which yield an effective ferromag- 
netic model. In order to make contact with previous work, 
we impose the following constraints : t a = 2t m , e a = 0, 
U arn = 2U aa and U aa = 2 — thus the energy scale is set 
by the on-site interaction. The parameters of the effective 
Ising Hamiltonian Eq. ^ are then found through the strong- 



coupling expansion 1271 



h = -2 + 6t*+e m + Ae m , (5) 
T = 2gV2. (6) 

We choose a representative value of £ a =5.04xl0 -2 
(J = 0.01), which avoids a slow down in the dynamics at 
lower values of t a . Eq. ^ predicts that the effective longi- 
tudinal field h is canceled by tuning e m = 1.9848. However, 
even after this, we observe in our simulations the effects of 
the presence of an effective longitudinal field, which arises 
from contributions in perturbation theory in small t a beyond 
the second-order expansion considered in Ref. 1211 . Rather 
than extending this calculation to higher order, we add an ex- 
tra perturbation, Ae m , to the molecular potential in Eq. ^ 
to compensate this "stray" longitudinal field. In general, this 




FIG. 1. (Color online) (a) Mapping between the pairing Bose-Hubbard 
model Eq. {TJ in the second Mott lobe and the Ising chain Eq. |2j. Two 
atoms (one molecule) per site map(s) to effective spin-down (spin-up). A lo- 
cal excitation converts one molecule into two atoms, which then propagate 
via Feshbach coupling and boson hopping (transverse field F and Ising in- 
teraction J, respectively), (b) Phase diagram of the full model Hamiltonian 
Eq. |TJ for t a = 5.04 x 10~ 2 (J = 0.01) and constraints detailed in the 
text, measured by the relative boson density, (c) Dynamical structure factor 
for 0=1.06 x 10~ 3 (r=0.3J) and Ae m =1.5 x 10~ 4 . The continuum of 
excitations reproduces exactly that of the Ising chain (inset), (d) Dynamical 
structure factor for g = 3.54 x 10 -3 (r = J), showing the same dispersion 
as the Ising chain in the paramagnetic state (inset). 



3 



term depends on the values of t a and g chosen. Probably, this 
is also a more realistic way of achieving the same result in an 
experimental setting, since it just involves further fine-tuning 
of e m , instead of balancing a version of Eq. ([5]) with more, 
higher-order, terms. 

The resulting phase diagram, where the Feshbach cou- 
pling g is the only free parameter and the order parameter is 
An=2S z , is shown in Fig.[TJb) for Ae m = 0. For low val- 
ues of g an ordered phase based on the product state with one 
molecule localised per site, |ff), is found. The stray longitudi- 
nal field described above biases the system towards this state, 
rather than the "atomic" ordered state based on two atoms lo- 
calised per site, All of the physical behaviour here de- 
scribed still holds in the case where the on-site potentials have 
been tuned to favour the atomic ground state. By increasing 
the value of g the system goes through a crossover into a disor- 
dered phase. This crossover is revealed by the non-diverging 
peak of the correlation length near g = 2 x 10 -3 (T « J/2), 
indicating that the Ising critical point is nearby. 

The excitation spectrum at low energies is revealed by the 
transverse dynamical structure factor function S y (k,uj), in- 
troduced above. The energy scales observed here are set by 
the effective Ising parameters and are hence rather low with 
uo « 0.01. This is well below the Mott gap U aa , above which 
single-particle excitations appear. We start with the "molecu- 
lar" ground state where the S y excitation in Eq. ^ dissociates 
one molecule into two atoms at site i, i.e. flips HJJ-). Evo- 
lution in time creates a domain of atoms, which is a bound 
state of two domain walls propagating in opposite directions 
along the chain. In the absence of a longitudinal field, the 
domain walls are deconfined as the two ground states are per- 
fectly degenerate. This is achieved by setting #=1.06 x 10 ~ 3 
(T=0.3 J) and fine tuning the on-site potential with a further 
term Ae m =1.5 x 10 -4 . The dynamical structure factor func- 
tion, shown in Fig.[TJc), displays a broad continuum of exci- 
tations around k = 0, sharpening up closer to the Brillouin 
zone edge at k = tt. This reproduces perfectly the behaviour 
of the same quantity in the transverse-field Ising chain, see 
e.g. Ref. 1321 , as shown in the inset, where the continuum of 
excitations arises from the incoherent movement of domain- 
wall pairs. A finite longitudinal field confines the domain 
walls, destroying the continuum of excitations, except for res- 
onances at specific momenta and energies, which can be seen 
as massive "meson" bound states l37l . An effective longitu- 
dinal field has the same effect on the molecular ground state 
under study, breaking up the continuum observed in Fig. [TJc) 
(not shown). 

Increasing g leads to the disordered region. The dynamical 
structure factor function is shown in Fig.[TJd) for a representa- 
tive value g=3.54 x 10 ~ 3 (T=J). A well-defined excitation 
with a quadratic dispersion relation is visible, in perfect cor- 
respondence with the spin-flip quasi-particle of the disordered 
Ising chain, whose structure factor is shown in the inset. 

The most impressive feature of the Ising model is the E$ 
symmetry 013, which emerges when the model is tuned to 
the critical point T=J/2 with a longitudinal field \h\ <C |T| 



applied. This is described by the perturbation of a c= 1/2 con- 
formal field theory, which is still integrable and manifests it- 
self in exactly eight massive bound states. The ratios between 
these masses, cf . Table [I] have been predicted analytically in 
Ref. |3 and observed in subsequent numerical studies, see e.g. 
Refs. I3ll32l. The mass ratio between the two lightest parti- 
cles is the golden ratio 1.618(...) and has been experimentally 
observed in CoM^Oe 0. However, heavier excitations are 
difficult to measure there, because their intensity is very close 
to that of the nearby continuum. It is not possible to mod- 
ify this by tuning the longitudinal field, which is fixed by the 
inter-chain coupling. 

m2/mi 7713/mi 77I4/77I1 ms/mi mfj/ 771 1 7777/7771 7778/7771 

1.618 1.989 2.405 2.956 3.218 3.891 4.783 

TABLE I. Analytically predicted mass ratios from Ref. |5 |. These 
numbers arise as eigenvalues of a matrix constructed from the roots 
of the Lie algebra Eg. 



The excitation spectrum when the system is tuned close 
to the perturbed Ising critical point, with g=1.80xl0 -3 
(r=0.51J)and Ae m =-5 x 10~ 4 (\h\ « |r|/10), is revealed 
by examining the dynamical structure factor function in 
Figs. |2ja,b). At least five different excitations are clearly 
identified above the continuum, which we can associate with 
the first four particles of E% theory and the bound- state pair 
mi + rri2- In Fig.[2jc) we present a sweep in Feshbach cou- 
pling near this point, while keeping the molecular potential 
perturbation fixed at Ae m =— 1 x 10 -3 for better convergence 
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FIG. 2. (Color online) Signatures of the emerging Eg symmetry in a Bose- 
Hubbard model, (a) Cut at k = of the dynamical structure factor function 
for g = 1.80 x 10~ 3 (T = J/2) and Ae m = -5 X 10~ 4 , showing several 
massive excitations, (b) Dynamical structure factor for the same parameters 
as (a) resolved in full momentum space, (c) Ratio of the masses of excitations 
to that of the lightest one, as a function of Feshbach coupling g, for fixed 
Ae m =— 1 x 10 -3 . Horizontal lines represent values analytically predicted 
from the Lie algebra E% . 
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of the results. The mass ratios calculated in this region are 
very close to the analytically predicted ones (horizontal lines), 
crossing them for a value of g = 1.83 x 10 -3 , which cor- 
responds to a transversal field T w 0.52 J, slightly larger 
than in the pure Ising chain. This small renormalisation of 
the critical value of field also probably arises from contribu- 
tions in higher order from perturbation theory. The differ- 
ent mass ratios increase roughly linearly with g. The heav- 
iest particles become progressively more difficult to observe 
with increasing g, since their spectral weight decreases, until 
the limit of Fig. [TJd) is eventually recovered for larger val- 
ues of g than the ones shown. The observation of this highly 
non-trivial sequence in the energy spectrum is a clear evi- 
dence for the emergence of E$ symmetry in a Bose-Hubbard 
model. The values obtained for the ratios do not change sig- 
nificantly when the molecular potential perturbation is dou- 
bled toAe m = —2 x 10 -3 , staying within error bars of the 
data points in Fig. [2jc). The stability of the results for such 
a relatively large ratio of fields |ft|/|r| « 1/2, has also been 
observed in the Ising chain |32l . 

Recent developments in Bragg spectroscopy applied to cold 
atomic systems, where two lasers with different frequencies 
shine on the sample at different angles, allow the study of 
the full excitation spectrum, resolved in momentum and en- 
ergy, even in the presence of an optical lattice (3U [39J. A 
complementary, and simpler, scheme of observing some of 
the behaviour described above in a cold-atoms setting is also 
desirable. Motivated by recent experimental [40 1 and theo- 
retical work EH [42), we now look at the propagation of a 
single excitation in space and time. We prepare the system 
in the molecular ground state for g=1.06 x 10 -3 , apply the 
S y dissociation excitation at a site io, and track the evolu- 
tion in time of the relative local density An at a site i lat- 
tice spacings apart. The effective longitudinal field has been 
tuned to h=0 by setting Ae m =1.5 x 10" 4 in Fig. [^a). At 
a qualitative level, An^ suddenly drops when a domain wall 
propagates through a site i, signalling the dissociation of the 
molecule (An>0) at that site into two atoms (An<0). Since 
there is no effective longitudinal field to confine the domain- 
walls, the An<0 domain persists in the long-time limit, and 
it grows with more sites being progressively flipped, lead- 
ing to the continuum of excitations in the dynamical struc- 
ture factor shown in Fig. [TJc). In Fig. [3jb), the behaviour 
in the presence of an effective longitudinal field is shown, 
for Ae m =— 1 x 10 -3 . The domain of atoms still grows up 
to a few sites away from the point where the excitation was 
added. However, it quickly collapses, due to the induced 
confining potential inhibiting the propagation of the domain 
walls. The sites which were excited eventually return to a state 
with An~l in the long-time limit, which is different from the 
original molecular ground state. 

Systems of hundreds of cold atoms confined to ID optical 
lattices have been extensively explored in the last decade 
[T8l . Recent experiments have investigated heteronuclear 
bosonic mixtures in ID lattices (43]|44]], which can be de- 
scribed by Bose-Hubbard models similar to Eq. ([T}, see 
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FIG. 3. (Color online) Time evolution of the relative density of bosons after 
a local quench dissociating two molecules into one atoms is performed on 
the molecular state for g=1.06 x 10 -3 . (a) By tuning Ae m =1.5 x 10 -4 , 
the effective longitudinal field is suppressed and the domain of atoms grows 
unbounded with time, (b) By tuning Ae m =— 1 x 10 -3 , a longitudinal field 
that thwarts the growth of the excited domain is induced. After a certain time 
the domain of atoms collapses, resulting in the confinement of the excitation. 

e.g. Ref. 1 45]. The two key requirements to observe the 
behaviour described here are the formation of a symmetry- 
breaking insulating phase and the presence of fluctuations able 
to destroy it. In order to ensure an insulating phase, the ratio 
of the hopping parameter to the intra-species potential, which 
depends on the lattice depth and the intra-species scattering 
length, should be t a /U aa < 0.1. Specific ratios between the 
different U s are overall not very important, leading only to the 
exchange value J. We have explicitly checked that our results 
do not change significantly for finite systems of comparable 
lengths L = 100 — 200. By the same reasoning, the presence 
of a trap potential should not affect much the results, as long 
as the central insulating domain with p T = 2 is large enough. 

The parameters that allow a direct tuning of the effective 
longitudinal and transverse fields of the Ising model Eq. ^ 
are, respectively, the Feshbach coupling g, which depends 
on the background scattering length and on the resonance 
width l23lL and the on-site potentials e a . Tuning the Feshbach 
coupling as in Fig. [3jc) allows a sweep near the critical point 
in order to find the optimal set of parameters corresponding to 
Eg symmetry. Biasing the on-site potentials gives a handle on 
the stability of the particles and their spectral weight relative 
to the nearby continuum. 

In conclusion, we have shown that the low-energy excita- 
tions of a ID Bose-Hubbard pairing model are well described 
by the dynamics of the quantum Ising chain. Our results indi- 
cate that the excitation spectrum of a ID system of cold atomic 
gases described by a model such as Eq. ([I]) can yield a very 
rich physical behaviour, from incoherent continuum of exci- 
tations to stable quasi-particles, culminating in the observa- 
tion of fingerprints of E% symmetry. Recent developments 
in manipulating heteronuclear bosonic mixtures in ID arrays, 
and in techniques directly accessing the low-energy excitation 
spectrum, open up the fascinating possibility of observing this 
behaviour in a cold atoms experiment. 
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